###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##
###*### TOOLBOX FUNCTIONS FOR:                          ###*###*###*###*###*###*###*###*##
###*### EQUILIBRIUM SIMULATING                          ###*###*###*###*###*###*###*###*##
###*### PARAMETER ESTIMATION & OUTCOME DISPLAY          ###*###*###*###*###*###*###*###*##
###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##

## Prints table results with friendly copy paste format for excel
print.x <- function(x) apply(x, 1, function(x){cat(x); cat("\n")})

## Vector of summary statistics for vector x
fsum <- function(x) c(mean(x), sd(x), min(x), max(x))

## Approximates the value of a constraint variable when out of bounds
bounds <- function(x, a = -1e+10, b = 1e+10, tol = 1e-6) {
  x <- min(max(a, x), b)
  return(x)
}

## Converts the mean and sd moments of a lognormal into its normal parameters
lnorpar <- function(mx, sx) {
  s <- sqrt(log(sx^2 / mx^2 + 1))
  m <- log(mx) - s^2/2
  #return(c('m' = m, 's' = s))
  return(c(m, s))
}

## Inverse of 'lnorpar'
lnorpar.inv <- function(m, s) {
  mx <- exp(m + s^2/2)
  sx <- sqrt((exp(s^2) - 1) * exp(2*m + s^2))
  return(c(mx, sx))
}

## Infer the implied correlation parameter of wage offer
lnorrho <- function(rhox, m1, m2, s1, s2) {
  mx1  <- exp(m1 + s1^2/2)
  mx2  <- exp(m2 + s2^2/2)
  cx12 <- exp(m1 + m2)*exp((s2^2 + (rhox + s1)^2)/2) - mx1*mx2
  rho  <- log(cx12/(mx1 * mx2) + 1) / (s1 * s2)
  return(as.numeric(rho))
}

## Inverse of 'lnorrho'
lnorrho.inv <- function(rho, m1, m2, s1, s2) {
  s12  <- rho * s1 * s2
  mx1  <- exp(m1 + s1^2/2)
  mx2  <- exp(m2 + s2^2/2)
  cx12 <- exp(m1 + m2 + (s1^2 + s2^2)/2) * (exp(s12) - 1)
  rhox <- sqrt(log((cx12 + mx1*mx2)/exp(m1 + m2))*2 - s2^2) - s1
  return(as.numeric(rhox))
}

## Transforms the parameter vector to apply constraints
conv <- function(parameters) {
  
  ## Child quality
  gm  <- parameters["gm"]  - parameters["gf"]
  gmh <- parameters["gmh"] - parameters["gm"]
  gmc <- parameters["gmc"] - parameters["gmh"]
  gR  <- parameters["gR"]  - parameters["gf"]
  gN  <- parameters["gN"]  - parameters["gR"]
  gL  <- parameters["gL"]  - parameters["gN"]
  gH  <- parameters["gH"]  - parameters["gL"]
  
  ## Lognormal parameters
  y1 <- lnorpar.inv(parameters['mwm'],  parameters['swm'])
  y2 <- lnorpar.inv(parameters['mwmh'], parameters['swmh'])
  y3 <- lnorpar.inv(parameters['mwmc'], parameters['swmc'])
  z1 <- lnorpar.inv(parameters['mcL'],  parameters['scL'])
  z2 <- lnorpar.inv(parameters['mcH'],  parameters['scH'])
  
  ## Wages offer
  mwmh <- y2["mwmh"] - y1["mwm"]
  mwmc <- y3["mwmc"] - y2['mwmh']
  rho  <- lnorrho.inv(parameters["rho"],  parameters['mwmh'], mwf, parameters['swmh'], swf)
  rhoc <- lnorrho.inv(parameters["rhoc"], parameters['mwmc'], mwf, parameters['swmc'], swf)
  
  ## Childcare costs
  mcH <- z2["mcH"] - z1["mcL"]
  
  pars <- c(parameters[c('ma1', 'ma2', 'ma3', 'sa1', 'sa2', 'sa3', 
                         'ma1s', 'ma2s', 'sa1s', 'sa2s', 'mu', 'tcc',
                         'g0', 'g0s', 'gf')], gm, gmh, gmc, parameters['gsm'],
            gL, gH, gN, gR, parameters['re'], y1, y2, y3, 
            'rho' = rho, 'rhoc' = rhoc, z1, z2)
  pars['mwmh'] <- mwmh
  pars['mwmc'] <- mwmc
  pars['mcH']  <- mcH
  
  return(pars)
  
}

## Reverts 'conv1'
conv.inv <- function(parameters, local = F) {
  
  parameters <- constraints(parameters)
  
  ## Household utility
  ma1  <- parameters["ma1"]
  ma2  <- parameters["ma2"]
  ma3  <- parameters["ma3"]
  sa1  <- parameters["sa1"]
  sa2  <- parameters["sa2"]
  sa3  <- parameters["sa3"]
  ma1s <- parameters["ma1s"]
  ma2s <- parameters["ma2s"]
  sa1s <- parameters["sa1s"]
  sa2s <- parameters["sa2s"]
  mu   <- parameters["mu"]
  tcc  <- parameters["tcc"]
  
  ## Child quality
  g0  <- parameters["g0"]
  g0s <- parameters["g0s"]
  gf  <- parameters["gf"]
  gm  <- parameters["gm"]  + gf
  gmh <- parameters["gmh"] + gm
  gmc <- parameters["gmc"] + gmh
  gsm <- parameters["gsm"]
  gR  <- parameters["gR"]  + gf
  gN  <- parameters["gN"]  + gR
  gL  <- parameters["gL"]  + gN
  gH  <- parameters["gH"]  + gL
  re  <- parameters["re"]
  
  ## Wages offer
  rho  <- parameters["rho"]
  rhoc <- parameters["rhoc"]
  mwm  <- parameters["mwm"]
  swm  <- parameters["swm"]
  mwmh <- parameters["mwmh"] + mwm
  swmh <- parameters["swmh"]
  mwmc <- parameters["mwmc"] + mwmh
  swmc <- parameters["swmc"]
  
  ## Childcare costs
  mcL <- parameters["mcL"]
  mcH <- parameters["mcH"] + mcL
  scL <- parameters["scL"]
  scH <- parameters["scH"]
  
  if(local) {
    mwm  <- mwm  * parameters['dwm']
    mwmh <- mwmh * parameters['dwm']
    mwmc <- mwmc * parameters['dwm']
    swm  <- swm  * parameters['dws']
    swmh <- swmh * parameters['dws']
    swmc <- swmc * parameters['dws']
    rho  <- rho  * parameters['drho'] 
    rhoc <- rhoc * parameters['drho'] 
    mcL  <- mcL  * parameters['dcm']
    mcH  <- mcH  * parameters['dcm']
    scL  <- scL  * parameters['dcs']
    scH  <- scH  * parameters['dcs']
  }
  
  ## Wage offer parameters
  y1   <- lnorpar(mwm,  swm)
  y2   <- lnorpar(mwmh, swmh)
  y3   <- lnorpar(mwmc, swmc)
  z1   <- lnorpar(mcL, scL)
  z2   <- lnorpar(mcH, scH)
  rho  <- lnorrho(rho, y2[1], mwf, y2[2], swf)
  rhoc <- lnorrho(rhoc, y3[1], mwf, y3[2], swf)
  
  pars <- c(ma1, ma2, ma3, sa1, sa2, sa3, ma1s, ma2s, sa1s, sa2s, mu, tcc,
            g0, g0s, gf, gm, gmh, gmc, gsm, gL, gH, gN, gR, re, 
            y1, y2, y3, 'rho' = rho, 'rhoc' = rhoc, 
            z1, z2)
  
  return(pars)
  
}

## Verifies and applies the predetermined constraints to the parameter vector
constraints <- function(parameters) {
  
  ## Non-negative parameters
  par  <- parameters
  ptve <- c("gf", "gm", "gmh", "gmc", "gL", "gH",
            "mcL", "mcH", "scL", "scH", "sa1", "sa2", "sa3", "mu",
            "mwm", "mwmh", "mwmc", "swm", "swmh", "swmc", "tcc")
  parameters[ptve]  <- abs(par[ptve])
  parameters['gsm'] <- -abs(par['gsm'])
  
  ## Boundary restrictions
  par                <- parameters
  parameters["rho"]  <- bounds(par["rho"],  a = -1,    b = 1)
  parameters["rhoc"] <- bounds(par["rhoc"], a = -1,    b = 1)
  parameters["mu"]   <- bounds(par["mu"],   a = 0.001, b = 1)
  parameters["re"]   <- bounds(par["re"],   a = -10,   b = 1)
  parameters["tcc"]  <- bounds(par["tcc"],  a = 0,     b = 3)
  # Preferences
  parameters["ma1"]  <- bounds(par["ma1"],  a = -1.100, b = 3.600)
  parameters["ma2"]  <- bounds(par["ma2"],  a = -1.100, b = 3.600)
  parameters["ma3"]  <- bounds(par["ma3"],  a = -1.100, b = 3.600)
  parameters["sa1"]  <- bounds(par["sa1"],  a = 0.001,  b = 3)
  parameters["sa2"]  <- bounds(par["sa2"],  a = 0.001,  b = 3)
  parameters["sa3"]  <- bounds(par["sa3"],  a = 0.001,  b = 3)
  parameters["ma1s"] <- bounds(par["ma1s"], a = -1.101, b = 3.600)
  parameters["ma2s"] <- bounds(par["ma2s"], a = -1.101, b = 3.600)
  parameters["sa1s"] <- bounds(par["sa1s"], a = 0.001,  b = 3)
  parameters["sa2s"] <- bounds(par["sa2s"], a = 0.001,  b = 3)
  # Quality prod.
  parameters["gf"]   <- bounds(par["gf"],  a = 0.01, b = 0.5)
  parameters["gm"]   <- bounds(par["gm"],  a = 0.01, b = 0.2)
  parameters["gmh"]  <- bounds(par["gmh"], a = 0.01, b = 0.2)
  parameters["gmc"]  <- bounds(par["gmc"], a = 0.01, b = 0.2)
  parameters["gsm"]  <- bounds(par["gsm"], a = -(parameters['gm'] + parameters['gf'])*0.95, b = 0)
  parameters["gL"]   <- bounds(par["gL"],  a = 0.01, b = 0.3)
  parameters["gH"]   <- bounds(par["gH"],  a = 0.03, b = 0.3)
  parameters["gR"]   <- bounds(par["gR"],  a = -parameters["gf"]+0.01, b = 0.2)
  # Costs
  parameters["scL"]  <- bounds(par["scL"], a = 0.0001, b = 5)
  parameters["scH"]  <- bounds(par["scH"], a = 0.0001, b = 5)
  # Wages
  parameters["mwm"]  <- bounds(par["mwm"],  a = 6.15,  b = 10)
  parameters["swm"]  <- bounds(par["swm"],  a = 0, b = 7)
  parameters["swmh"] <- bounds(par["swmh"], a = 0, b = 7)
  parameters["swmc"] <- bounds(par["swmc"], a = 0, b = 14)
  
  return(parameters)
  
}

## Criterion function to be minimized for parameter estimation
crit_fun <- function(parameters, sample_moments, weights = 1) {
  sim_moments <- simul_market(parameters)$smoms_fin
  dist        <- (sim_moments - sample_moments)*weights
  i           <- which(!is.na(sample_moments) & is.na(sim_moments))
  dist[i]     <- 1000
  dist        <- mean(dist^2, na.rm = TRUE)
  #print(dist)
  return(dist)
}

## Allows to transform crit_fun to only depend of a subset of parameters
crit_funn <- function(px, parameters, nampar, sample_moments, weights = 1) {
  print(round(px, 6))
  parameters[nampar] <- px
  dist <- crit_fun(parameters, sample_moments, weights)
  print(dist)
  return(dist)
}

## Market structure
gstate <- function(k, gL, gH) {
  # Defines the k-th vector of center quality available in the market
  # out of all the k = 1,...,3^NL possible combinations
  gSpace <- c(0, gL, gH)
  nx     <- length(gSpace)
  x <- sapply(1:NL, function(i) {
    ik <- ceiling((k - nx^(i)*floor(k/nx^(i)))/(nx^(i-1)))
    if(ik == 0) ik <- nx
    gSpace[ik]
  })
  return(x)
}

## Child quality vector
ChildQ <- function(q0, ta, tb, gd, g0, gf, gm, gmh, gmc, gsm, re, ed) {
  gmi <- gm*(ed == 1) + gmh*(ed == 2) + gmc*(ed == 3) + gsm*sm
  q   <- q0 + g0 + (gf*ta^re + gmi*tb^re + gd*(8^re))^(1/re)
  return(q)
}

## Child quality scalar
childQ <- function(q0, ta, tb, gd, g0, gf, gmi, re) {
  q   <- q0 + g0 + (gf*ta^re + gmi*tb^re + gd*(8^re))^(1/re)
  return(q)
}

## Utility as fn of ta
uta <- function(ta, ha, hb, gd, q0, a2, a3, a4, g0, gf, gmi, re, gL, gH, gN, tcc) {
  # Version of the utility when all variables but ta are fixed.
  # Used to solve the problem of dad's time allocation, given a mom's work choice
  # and a child care arrangement.
  cc <- hb == 4 & gd == 0
  cd <- hb == 8 & gd != gN
  ta <- max(0.01 + 4*cc + tcc*cd, min(7.99, ta))
  tb <- 16 - ta - 8*(gd > 0)
  q  <- childQ(q0, ta, tb, gd, g0, gf, gmi, re)
  la <- 16 - ha - ta
  lb <- 16 - hb - tb - tcc*(gd != gN & gd > 0)
  un <- suppressWarnings(as.numeric(a2*log(q) + a3*log(la) + a4*log(lb)))
  if(is.nan(un) | is.infinite(un)) un <- 100000
  return(-un)
}

## Organizes ta as matrix
FindTstar <- function(tstar, g, hhopts) {
  # Finds the optimal ta for every given mom's labor choice and child care arrangement
  ta <- matrix(0, nrow(tstar), 3*(NL + 2) + 2)
  ta[, 3*(NL + 2) + 1:2] <- tstar[, c("8,4,0", "8,0,0")]
  # print(sort(g))
  # print(sort(unique(hhopts$gd)))
  for(i in which(g > 0)) ta[, 1:3 + 3*(i - 1)] <- tstar[, which(hhopts$gd == g[i])]
  return(ta)
}

## Household's utility function
HHUtility <- function(p, gmat, q, la, lb, inco, c_disc, smdisc, rel, a1, a2, a3, a4) {
  # Child care price arrangement
  pmat <- cbind(matrix(rep(p, each = 3), nrow(gmat), 3*(NL + 2), byrow = TRUE), 0, 0)
  pmat[-rel, -2:0 + 3*(NL + 2)] <- NA
  pmat <- pmat*(1 - smdisc)
  pmat <- pmat - c_disc
  pmat[pmat < 0] <- 0
  # Consumption
  cons <- inco - 8*(gmat > 0)*pmat
  cons[cons <= 0] <- NA
  # Mean utility (before the preferences shock)
  u  <- suppressWarnings(a1*log(cons) + a2*log(q) + a3*log(la) + a4*log(lb))
  return(list(u = u, cons = cons))
}

## Market shares demand
Shares <- function(u, g, mu, rel) {
  # Defines the probabilities that a given child care arrangement is hired
  # given the model parameters and a market supply structure.
  open <- 1*(g > 0)
  open <- cbind(matrix(rep(open, each = 3), nrow(u), 3*(NL + 2), byrow = TRUE), 1, 1)
  open[-rel, -2:0 + 3*(NL + 2)] <- 0
  open[sm == 1, ncol(open)-1] <- 0
  # Define probabilities by household type
  num  <- exp(u/mu)*open
  num[is.na(num)]   <- 0
  num[num > 1e+200] <- 1e+201
  den   <- matrix(rowSums(num), nrow(u), 3*(NL + 2) + 2)
  probs <- num/den
  probs[num == 0] <- 0
  m     <- seq(1, 3*(NL + 2), by = 3)
  # Define market shares (expected demand) as the average by alternative
  s     <- colMeans(probs[, m] + probs[, m + 1] + probs[, m + 2])
  return(list(probs = probs, s = s))
}

## Profits of open firms
Profits <- function(p, g, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                    a1, a2, a3, a4, mu) {
  # Utility and consumption
  U <- HHUtility(p, gmat, q, la, lb, inco, c_disc, smdisc, rel, a1, a2, a3, a4)
  cons <- U$cons
  u  <- U$u
  # Shares and demand
  s  <- Shares(u, g, mu, rel)
  N  <- s$s*nrow(u)*M
  # Profits
  pi <- (p - mc)*N*8 - Fi*ceiling(N/40)
  pi <- pi*(g > 0)
  return(list(out = data.frame(g = g, p = p, pi = pi, mc = mc, s = s$s),
              u = u, cons = cons, probs = s$probs))
}

## Predicted markups
Markups <- function(p, g, gmat, q, la, lb, inco, c_disc, smdisc, 
                    rel, a1, a2, a3, a4, mu) {
  # Computes the firms' optimal predicted markups from FOC
  # Utility
  HH     <- HHUtility(p, gmat, q, la, lb, inco, c_disc, smdisc, rel, a1, a2, a3, a4)
  cons   <- HH$cons
  u      <- HH$u
  # Probabilities and shares
  shares <- Shares(u, g, mu, rel)
  s      <- shares$s
  sj     <- shares$probs
  # Derivatives
  m      <- seq(1, 3*(NL + 2), by = 3)
  den    <- (8 * sj * (1 - sj) * a1/mu) / cons
  den    <- colMeans(den[, m] + den[, m + 1] + den[, m + 2], na.rm = TRUE)
  # Markup
  mkup         <- s / den
  mkup[NL + 2] <- 0
  mkup[s == 0] <- 0
  return(mkup)
}

## Fixed point to find static eq using markups
FP  <- function(p, g, gmat, mc, q, la, lb, inco, c_disc, smdisc, rel, 
                a1, a2, a3, a4, mu) {
  # Compares the input price, p, with the optimal price implied by FOC eq. (16)
  # p_new, to find the market equilibrium
  # Markups
  mkup  <- Markups(p, g, gmat, q, la, lb, inco, c_disc, smdisc, rel, 
                   a1, a2, a3, a4, mu)
  p_new <- mc + mkup
  ## Shares
  HH     <- HHUtility(p_new, gmat, q, la, lb, inco, c_disc, smdisc, rel, a1, a2, a3, a4)
  shares <- Shares(HH$u, g, mu, rel)
  s      <- shares$s
  # Evaluate difference
  dist  <- sqrt(mean((p - p_new)^2, na.rm = TRUE))
  if(is.na(dist)) dist <- 10000
  return(list(dist = dist, p_new = p_new))
}

## Static equilibrium
Profeq2 <- function(gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                    a1, a2, a3, a4, mu) {
  # Algorithm to apply the fixed point function, FP and find the equilibrium 
  # price and demand
  p     <- p0   <- (mc + c(rep(2, NL), rep(0, 2)))*(gs > 0)
  iter2 <- iter <- 0
  aw    <- 0.2
  crit0 <- 1000
  # Iterate over p
  repeat {
    iter <- iter + 1
    ans  <- FP(p0, gs, gmat, mc, q, la, lb, inco, c_disc, smdisc, rel, 
               a1, a2, a3, a4, mu)
    crit <- as.numeric(ans$dist)
    p0   <- ans$p_new*aw + p0*(1 - aw)
    p0[is.na(p0)] <- p[is.na(p0)]
    if(crit < crit0) {
      iter2 <- 0
      crit0 <- crit
      p     <- ans$p_new
      #print(c(crit, aw))
    } 
    iter2 <- iter2 + 1
    if(iter2 > 10) aw <- min(max(0.05, aw/crit), 1) 
    if(crit < 0.01 | iter > 200) break
  }
  # Equilibrium profits
  ph  <- p
  out <- Profits(ph, gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                 a1, a2, a3, a4, mu)
  # Equilibrium output
  return(out)
}

## Fixed point to find: p = MeC
ZeroProf_PP  <- function(mk, gs, gmat, mc, q, la, lb, inco, c_disc, smdisc, rel, 
                         a1, a2, a3, a4, mu, Fi) {
  # For the public provision experiment: computes profits for a given markup (mk)
  # Non-negative function that should be minimized to find the closest approximation
  # to zero profits.
  mk <- c(abs(mk), numeric(NL - length(mk) + 2))
  p  <- mc + mk
  x  <- Profits(p, gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, 
                rel, a1, a2, a3, a4, mu)
  pi  <- x$out$pi[1:length(mk)]
  pi  <- pi * (1 + 1*(pi < 0))
  out <- sqrt(sum(pi^2))
  # print(round(mk, 2))
  # print(round(sum(out), 2))
  return(out)
}

## Static equilibrium: Public Provision
Profeq_PP <- function(gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                      a1, a2, a3, a4, mu) {
  # For the public provision experiment: finds the price that results in 
  # zero profits by minimizing ZeroProf_PP for the given set of open center providers.
  mk0 <- numeric(sum(gs[1:NL] > 0))
  ans <- optim(mk0, fn = ZeroProf_PP, gs = gs, gmat = gmat, mc = mc, q = q, 
               la = la, lb = lb, inco = inco, c_disc = c_disc, smdisc = smdisc, 
               rel = rel, a1 = a1, a2 = a2, a3 = a3, a4 = a4, mu = mu, Fi = Fi, 
               control = list(maxit = 500))
  mk  <- ans$par
  mk  <- c(abs(mk), numeric(NL - length(mk) + 2))
  ph  <- mc + mk
  # Equilibrium profits
  out <- Profits(ph, gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                 a1, a2, a3, a4, mu)
  # sum(out$out$pi)
  # Equilibrium output
  return(out)
}

## Competition stage equilibrium outcomes
master <- function(gx, Cg, Cn, Cr, Fixed, gL, gH, gN, gR, ha, hb, 
                   q0, g0, gf, gm, gmh, gmc, gsm, re, tstar, hhopts, 
                   ed, inco, c_disc, smdisc, hqv, 
                   rel, a1, a2, a3, a4, mu, tcc) {
  # Computes all the outcomes to be used of the competition stage of the supply-side
  # game, given the parameters, households' time allocation arrangements and a 
  # supply market structure of child care quality availability.
  g    <- as.numeric(gx) # Firm quality arrangement
  gs   <- c(g, gN, gR)
  gmat <- matrix(c(rep(gs, each = 3), 0, 0), nrow(hb), 3*(NL + 2) + 2, byrow = TRUE)
  gmat[-rel, 1:3 + 3*(NL + 2) - 3]  <- 0
  # Quality voucher
  if(hqv) c_disc[, gmat[1, ] != gH] <- 0
  # Marginal costs
  mc <- c((g == gL)*Cg[, 1] + (g == gH)*Cg[, 2], Cn, Cr)
  Fi <- Fixed[1]*(gs == gL) + Fixed[2]*(gs == gH)
  # Time allocation
  ta <- FindTstar(tstar, gs, hhopts)
  td <- cbind(matrix(8, J, 3*(NL + 2)), 0, 0)
  tb <- 16 - ta - td
  la <- To - ha - ta 
  lb <- To - hb - tb - tcc*(gmat != gN & gmat > 0)
  lb[sm == 1, ncol(lb)-1] <- NA
  # Child quality for each scenario in gmat
  q   <- ChildQ(q0, ta, tb, gmat, g0, gf, gm, gmh, gmc, gsm, re, ed)
  # Compute 2nd-stage equilibrium (conditional on g)
  out <- Profeq2(gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, rel, 
                 a1, a2, a3, a4, mu)
  return(out)
}

## Competition stage equilibrium outcomes: PP
master_PP <- function(gx, Cg, Cn, Cr, Fixed, gL, gH, gN, gR, ha, hb, 
                      q0, g0, gf, gm, gmh, gmc, gsm, re, tstar, hhopts, 
                      ed, inco, c_disc, smdisc, hqv, 
                      rel, a1, a2, a3, a4, mu, tcc) {
  # For the public provision experiment: Computes all the outcomes to be used 
  # from the competition stage of the supply-side game, given the parameters, 
  # households' time allocation arrangements and a  supply market structure of 
  # child care quality availability.
  g    <- as.numeric(gx) # Firm quality arrangement
  gs   <- c(g, gN, gR)
  gmat <- matrix(c(rep(gs, each = 3), 0, 0), nrow(hb), 3*(NL + 2) + 2, byrow = TRUE)
  gmat[-rel, 1:3 + 3*(NL + 2) - 3]  <- 0
  # Quality voucher
  if(hqv) c_disc[, gmat[1, ] != gH] <- 0
  # Marginal costs
  mc <- c((g == gL)*Cg[, 1] + (g == gH)*Cg[, 2], Cn, Cr)
  Fi <- Fixed[1]*(gs == gL) + Fixed[2]*(gs == gH)
  # Time allocation
  ta <- FindTstar(tstar, gs, hhopts)
  td <- cbind(matrix(8, J, 3*(NL + 2)), 0, 0)
  tb <- 16 - ta - td
  la <- To - ha - ta 
  lb <- To - hb - tb - tcc*(gmat != gN & gmat > 0)
  lb[sm == 1, ncol(lb)-1] <- NA
  # Child quality for each scenario in gmat
  q   <- ChildQ(q0, ta, tb, gmat, g0, gf, gm, gmh, gmc, gsm, re, ed)
  # Compute 2nd-stage equilibrium (conditional on g)
  out <- Profeq_PP(gs, gmat, mc, Fi, q, la, lb, inco, c_disc, smdisc, 
                   rel, a1, a2, a3, a4, mu)
  return(out)
}

## Display the simultation results
cf_table <- function(m) {
  # Function to display the equilibrium results out of the 
  # simul_market and simul_market_PP functions
  
  # Mental socres
  mt   <- weighted.mean(m$oq, m$share_new)
  i    <- which(sm == 0 & ed == 3)
  mm3  <- weighted.mean(m$oq[i, ], m$share_new[i, ])
  i    <- which(sm == 0 & ed == 2)
  mm2  <- weighted.mean(m$oq[i, ], m$share_new[i, ])
  i    <- which(sm == 0 & ed == 1)
  mm1  <- weighted.mean(m$oq[i, ], m$share_new[i, ])
  i    <- which(sm == 1)
  msm  <- weighted.mean(m$oq[i, ], m$share_new[i, ])
  # Non-Relative care
  k    <- which(m$gmat[1, ] == m$pars['gN'])
  pNt  <- sum(m$share_new[, k])/sum(m$share_new)
  i    <- which(sm == 0 & ed == 3)
  pNm3 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 2)
  pNm2 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 1)
  pNm1 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 1)
  pNsm <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  # Relative care
  k    <- which(m$gmat[1, ] == m$pars['gR'])
  pRt  <- sum(m$share_new[, k])/sum(m$share_new)
  i    <- which(sm == 0 & ed == 3)
  pRm3 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 2)
  pRm2 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 1)
  pRm1 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 1)
  pRsm <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  # HQ center care
  k     <- which(m$gmat[1, ] == m$pars['gH'])
  pHQt  <- sum(m$share_new[, k])/sum(m$share_new)
  i     <- which(sm == 0 & ed == 3)
  pHQm3 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 0 & ed == 2)
  pHQm2 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 0 & ed == 1)
  pHQm1 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 1)
  pHQsm <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  # LQ center care
  k     <- which(m$gmat[1, ] == m$pars['gL'])
  pLQt  <- sum(m$share_new[, k])/sum(m$share_new)
  i     <- which(sm == 0 & ed == 3)
  pLQm3 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 0 & ed == 2)
  pLQm2 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 0 & ed == 1)
  pLQm1 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i     <- which(sm == 1)
  pLQsm <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  # Center
  cER  <- sum(m$SPEQ$g[1:NL] > 0)/NL*100
  cERH <- sum(m$SPEQ$g[1:NL] == m$pars['gH'])/NL*100
  cERL <- sum(m$SPEQ$g[1:NL] == m$pars['gL'])/NL*100
  avsh <- mean(m$SPEQ$s[1:NL][m$SPEQ$g[1:NL] > 0])*100
  cpf  <- mean(ceiling(m$SPEQ$s[m$SPEQ$g %in% m$pars[c('gH', 'gL')]]*J*M/40))
  cpfH <- mean(ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gH']]*J*M/40))
  cpfL <- mean(ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gL']]*J*M/40))
  cn   <- sum(ceiling(m$SPEQ$s[1:NL]*J*M/40))
  cnH  <- sum(ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gH']]*J*M/40))
  cnL  <- sum(ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gL']]*J*M/40))
  cAE  <- mean(m$SPEQ$s[m$SPEQ$g %in% m$pars[c('gH', 'gL')]]*J*M / ceiling(m$SPEQ$s[m$SPEQ$g %in% m$pars[c('gH', 'gL')]]*J*M/40))
  cAEH <- mean(m$SPEQ$s[m$SPEQ$g == m$pars['gH']]*J*M / ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gH']]*J*M/40))
  cAEL <- mean(m$SPEQ$s[m$SPEQ$g == m$pars['gL']]*J*M / ceiling(m$SPEQ$s[m$SPEQ$g == m$pars['gL']]*J*M/40))
  cAP  <- weighted.mean(m$SPEQ$p[1:NL], w = m$SPEQ$s[1:NL])
  cAPH <- weighted.mean(m$SPEQ$p[m$SPEQ$g == m$pars['gH']], w = m$SPEQ$s[m$SPEQ$g == m$pars['gH']])
  cAPL <- weighted.mean(m$SPEQ$p[m$SPEQ$g == m$pars['gL']], w = m$SPEQ$s[m$SPEQ$g == m$pars['gL']])
  # Work moms
  k    <- which(m$hb[1, ] > 0)
  lpt  <- sum(m$share_new[, k])/sum(m$share_new)
  i    <- which(sm == 0 & ed == 3)
  lpm3 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 2)
  lpm2 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 0 & ed == 1)
  lpm1 <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  i    <- which(sm == 1)
  lpsm <- sum(m$share_new[i, k])/sum(m$share_new[i, ])
  # Wages
  k   <- which(m$hb[1, ] > 0)
  wt  <- weighted.mean(m$wi[, 2], w = rowSums(m$share_new[, k]))
  i   <- which(sm == 0 & ed == 3)
  wm3 <- weighted.mean(m$wi[i, 2], w = rowSums(m$share_new[i, k]))
  i   <- which(sm == 0 & ed == 2)
  wm2 <- weighted.mean(m$wi[i, 2], w = rowSums(m$share_new[i, k]))
  i   <- which(sm == 0 & ed == 1)
  wm1 <- weighted.mean(m$wi[i, 2], w = rowSums(m$share_new[i, k]))
  i   <- which(sm == 1)
  wsm <- weighted.mean(m$wi[i, 2], w = rowSums(m$share_new[i, k]))
  # Output
  out <- c(mt, mm3, mm2, mm1, msm, 
           c(pNt, pNm3, pNm2, pNm1, pNsm, pRt, pRm3, pRm2, pRm1, pRsm, 
             pHQt, pHQm3, pHQm2, pHQm1, pHQsm, 
             pLQt, pLQm3, pLQm2, pLQm1, pLQsm)*100, 
           cER, cERH, cERL, avsh, cpf, cpfH, cpfL, cn, cnH, cnL, 
           cAE, cAEH, cAEL, cAP, cAPH, cAPL,
           c(lpt, lpm3, lpm2, lpm1, lpsm)*100, wt, wm3, wm2, wm1, wsm)
  return(out)
}
